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With the wide availability of affordable multiple-core parallel supercomputers, next generation 
numerical simulations of flow physics are being focused on unsteady computations for problems 
involving multiple time scales and multiple physics. These simulations require higher solution 
accuracy than most algorithms and computational fluid dynamics codes currently available. This 
paper focuses on the developmental effort for high-fidelity multi-dimensional, unstructured-mesh 
flow solvers using the space-time conservation element, solution element (CESE) framework. Two 
approaches have been investigated in this research in order to provide high-accuracy, cross-cutting 
numerical simulations for a variety of flow regimes: 1) time-accurate local time stepping and 2) high- 
order CESE method. The first approach utilizes consistent numerical formulations in the space-time 
flux integration to preserve temporal conservation across the cells with different marching time 
steps. Such approach relieves the stringent time step constraint associated with the smallest time 
step in the computational domain while preserving temporal accuracy for all the cells. For flows 
involving multiple scales, both numerical accuracy and efficiency can be significantly enhanced. The 
second approach extends the current CESE solver to higher-order accuracy. Unlike other existing 
explicit high-order methods for unstructured meshes, the CESE framework maintains a CFL 
condition of one for arbitrarily high-order formulations while retaining the same compact stencil as 
its second-order counterpart. For large-scale unsteady computations, this feature substantially 
enhances numerical efficiency. Numerical formulations and validations using benchmark problems 
are discussed in this paper along with realistic examples. 
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length of the cavity 

frequency of instability oscillations 

freestream Mach number 

Reynolds number based on length 

cylinder height for both protuberance and cavity 

depth of the cavity, or diameter of the cylindrical roughness element 

boundary-layer thickness 

time 

streamwise coordinate 
wall-normal coordinate 

spanwise coordinate in direction parallel to leading edge 
density 

streamwise velocity 
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V 


= wall-normal velocity 
w = spanwise velocity 

e = total energy, defined in the second section 

p = pressure 

R = gas constant 

T = temperature 

U = dependent solution vector, defined in the second section 

V = integration volume 

£2 = surface for the integration volume 

F = flux vectors for general conservation laws 

F x F y F z = flux vectors for compressible Navier-Stokes equations in three spatial directions 

S = source vector for general conservation laws 

(9, Vi = solution points and vertices within each solution element 

V! = stretched vertices to add numerical dissipation in the EBD approach 

Amin = global minimum spacing between two vertices 

Qmm Gmw = global lower and upper limit of the numerical dissipation parameter 

jc 0 , ,yo , zo = coordinates of solution point within the solution element 

y = ratio of the specific heat 

h = flux in the joint space-time domain 

s = surface normal vector in the joint space-time domain 

A = area of the space-time element interface 

d p = shock detection parameter 

s = small number 1 0" 20 

Qi = magnitude of derivative vectors 

Ci = sorted values of 0 t 

y/i = intermediate variables used to derive weighted average 

Wt = weighting for i-th derivative 

a = coefficient of the scalar wave equation 

0 = wave variable used to define solutions to the scalar wave equation 

A w = acoustic wave amplitude 

co = angular frequency of acoustic wave 

k x> k y , k z = wave numbers in jc, y, and z direction, respectively 

Subscripts 

x, y,z,t = derivatives in spatial and temporal directions 

Superscripts 

1 = number index for surrounding element and derivatives 


I. Introduction 

A fter decades of computational fluid dynamics (CFD) research, the primary development in recent 
years is focused on high-fidelity numerical computations for complex configurations. While 
confidence in predictive capabilities using CFD is relatively high in certain domains such as transonic 
aerodynamics using steady-state Reynolds Averaged Navier-Stokes (RANS) equations 1 , available codes 
and numerical algorithms to-date are still deficient in handling time accurate computations over complex 
configurations with multiple physical scales. For example, computations of unsteady wakes behind a large 
roughness element in supersonic or hypersonic boundary layers 2 * " 4 remain quite challenging. The 
spontaneous instability waves due to the presence of the roughness element could trigger unstable, high- 
frequency waves downstream that eventually could lead to laminar-to-turbulent transition. Several issues 
are associated with such simulations. First, to handle the geometry with a reasonable grid count, an 
unstructured mesh is favored, despite the fact that most matured and well-used CFD codes are still 
structured-mesh based. Second, temporal accuracy is a main concern. To simulate such highly unstable 
flows, numerical dissipation has to be sufficiently reduced in order not to hamper the generation and 
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development of instability waves. Third, computation of transitional flows or fully turbulent flows without 
any modeling requires resolution of multiple temporal and spatial scales. The challenging aspect is to 
preserve accuracy with large disparity of mesh sizes while maintaining computational efficiency. Last, but 
not least, is the requirement to handle interactions of fine scale waves with flow discontinuities. Shock 
capturing inevitably demands larger numerical dissipation that could in turn improperly suppress the wave 
development. Improving numerical algorithms and developing efficient parallel computer codes to meet 
the above requirements are no doubt one of the most important research directions for CFD development in 
the next decade. 

The importance of using unstructured meshes for complex configurations goes beyond doubts. 
Complex geometries post great challenges for structured grid generation. Laying out topological 
arrangement to cope with corners, small parts, and minute details appears to be quite time-consuming. 
Even with the success of the topology layout, connectivity requirements among blocks inevitably add 
unwanted clustering that in turn substantially increases the total grid count. Unstructured mesh generation 
on the other hand is free of these issues. An unstructured mesh also makes grid adaptation easier to 
implement because it does not have topological constraints. To maintain high temporal accuracy using 
unstructured meshes, several high-order numerical algorithms have been investigated in recent years. The 
discontinuous Galerkin (DG) method formulated by a combination of the weak form and discontinuous 
cell-interface integrals to allow jump solutions has been actively pursued. 5 " 10 The spectral difference 
method 11 " 12 also offers an alternative way to attain high-order accuracy for unstructured meshes. While 
most of these high-order methods offer different formulations to capture shocks and treat viscous terms, 
they follow a similar strategy to calculate cell-interface fluxes via an approximate ID Riemann solver. 
Special attention is often required to deal with multi-dimensional flows with shocks. 

The space-time conservation element, solution element (CESE) method introduced by Chang 13 in the 
mid 1990s offers a distinctly different numerical framework to solve general conservation laws via a 
consistent formulation over the space-time domain. Unlike conventional methods following the ID 
Riemann approximation concept, the CESE method is a genuine multi-dimensional scheme formulated by 
enforcing space-time flux conservation in the discretized domain. Integration of the conservation laws is 
carried out on staggered discretized volumes called conservation elements (CE), whose boundaries 
encompass only the smooth, well-defined regions defined by a solution element (SE). Consequently, no 
interface reconstruction or dimensional-splitting is needed. Flux conservation in space and time offers 
distinctly better numerical properties, both in terms of boundary condition treatment and temporal accuracy. 
It has been shown that this relatively newer numerical framework handles waves and shocks with accuracy 
comparable to existing high-order methods. 14 ' 15 The CESE method consists of a family of numerical 
schemes for the solutions of general conservation laws. These schemes form the foundation of the CESE 
numerical framework to be discussed in this paper. 

With the wide availability of parallel supercomputers, numerical simulations of flow physics will be 
focused on time accurate computations with multiple time scales and physics. This new demand stems from 
the fact that most outstanding flow problems involve unsteady governing equations along with a large 
disparity in time scales in different parts of the flowfield. While modeling such intricate physical 
phenomena remains important and in fact, is still the only feasible way for many cases, direct time-accurate 
computations play an increasingly important role. This paper is focused on developing unstructured mesh 
based numerical algorithms and computer software framework with high temporal accuracy for unsteady 
flow simulations. To this end, numerical algorithm studies are carried out in this research to provide high- 
fidelity, cross-cutting numerical simulations for a variety of flow regimes via two main approaches: 1) 
time-accurate local time stepping and 2) high-order CESE method. Both investigations are part of the 
continual development of the NASA Langley Research Center’s unstructured-mesh Navier-Stokes solver 
ez4d as part of the research effort under the Revolutionary Computational Aerosciences (RCA) subproject 
of the Fundamental Aeronautical Program (FAP). 

Time-accurate numerical computations can be done with either implicit or explicit methods. Explicit 
methods have a time step constraint with a CFL number around one as a direct outcome of the finite-size 
domain of dependence. Implicit methods, on the other hand, maintain numerical stability at a much larger 
CFL number. Depending on the time scale of the flow physics, implicit methods may or may not offer 
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higher computational efficiency. For low speed problems with a time scale much larger than the acoustic 
wave time scale, implicit methods are much more efficient. For supersonic flows, the instability waves and 
near wall structures have a time scale comparable to or even smaller than the acoustic wave scale. 
Accuracy requirements thus dictate a small time step around a CFL number of one. Many direct numerical 
simulations of high-speed flows are performed with an explicit time-marching algorithm. For these kinds 
of problems, an implicit method may still be preferred due to its robustness but not accuracy considerations. 
The CESE framework is formulated as an explicit method and the stable CFL limit is one, similar to other 
widely used explicit methods for unstructured meshes. For the CESE method, numerical dissipation 
increases as the CFL number decrease. 13 For a computational domain with multiple temporal-spatial scales, 
a constant time step results in a large disparity in CFL number. Regions with large grid size become less 
accurate if a constant time step is used throughout the domain. A time-accurate local time-stepping 
algorithm with a uniform CFL number while preserving time-accuracy could significantly improve both 
numerical accuracy and computational efficiency. Parallel to the development of local time-stepping, a 
high-order CESE method could also increase solution accuracy for time-accurate computations with 
multiple scales. Numerical algorithm advancement in these two areas is aimed at providing a more 
accurate flow solver for high-fidelity flow physics simulations. 

The time accurate local time-stepping implementation in the full 3D Navier-Stokes CESE solver 
follows the concept proposed by Chang et al. 16 for one-dimensional flows. The key idea is to derive the 
local time-stepping algorithm by enforcing flux conservation both in space and in time. Since the CESE 
method is built upon flux conservation over the entire space and time domain, the implementation of 
different time steps at different cells is relatively straightforward. A book keeping procedure similar to 
what was done by Yen 17 is adopted. The time steps across different cells are divided in a set of pre- 
determined discrete values. Depending on the time steps of the neighboring cells, the solutions are stored 
at multiple levels. These stored solutions are then used for the flux integration along the CE interfaces to 
make sure flux conservation is enforced at all times. For the high-order CESE development, the method 
proposed by Chang 18 is extended for 3D Navier-Stokes equations. In contrast to conventional high-order 
DG methods, the high-order CESE schemes proposed by Chang maintain a CFL limit of one regardless of 
the order of accuracy. In contrast, the stable CFL number limit decreases as the order of polynomials 
increases for explicit DG methods. For example, the stable CFL limit is only 0.1 for the p4 scheme. 19 A 
CFL limit of one for the high-order CESE method offers much relief from this stringent time step limitation 
as compared to other explicit methods. The implementation of Chang’s high-order CESE method by 
Bilyeu et al. 20 " 21 verifies the high CFL limit for 2D Euler equations. They also show that it is possible to 
construct higher than the 4th-order accurate CESE schemes for ID wave equations. Both time accurate 
local time-stepping schemes and the 4th-order accurate CESE framework are implemented in the ez4d code 
for large-scale 3D Navier-Stokes Computations. Numerical algorithms implemented are discussed in the 
next section, followed by results and discussion. Some concluding remarks are included at the end of this 
paper. 


II. Numerical Method 


Three-dimensional general conservation laws can be written as 


dU 

dt 


+ V'F = S 


( 1 ) 


where U is the dependent variable, F is the flux vector, and S is the source term. Both F and S are 
functions of dependent variables. For the special case of compressible Navier-Stokes equations, the 
dependent variables are U = (p, pu , pv, pw , e) and the flux vectors F = (F X9 F y , F z ) satisfy the following 
differential equation 


dU dF dF dF 

+ — - + — - + — L 

dt dx dy dz 


=0 


( 2 ) 
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where x, y, z, and t represent spatial coordinates and time, respectively. Flow variables p, u , v, w, and e 
represent density, three velocity components, and total energy (e = p /(/ — 1) + p(u 2 + V 2 + W 2 )/ 2 

where p denotes the pressure), respectively. Definitions of the total flux vectors (including an Euler part 
and a viscous part) F x , F y , and F z , can be found in standard text books and will not be included here. The 
source vector, S, contains all external forcing or other energy-related source terms; and it vanishes for 

compressible Navier-Stokes equations. To close the system, the perfect gas relation, p = pRT , with T 
representing the temperature and R the gas constant, is used in conjunction with eq. (1). 

Numerical solutions of eq. (1) can be obtained either by a differential or integral form. The integral 
form is preferred for fluid dynamics equations because of two main reasons. First, the integral form 
enforces flux conservation in the discretized space. Second, numerical issues associated with derivatives 
across flow discontinuities can be treated consistently only in the integral form. Most modem CFD solvers 
are formulated in the integral form but with distinct differences in the way eq. (1) is integrated. Most 
existing upwind or DG formulations treat the temporal derivative in eq. (1) as a separate identity and a 
finite-difference form is assumed for temporal derivatives. As a result, flux conservation is only enforced 
in the spatial domain. In the discretized space, this is equivalent to introducing a temporal source term in 
the conservation law. This source term depends on temporal discretization as well as the time step 
employed. While this in general works for solving steady state, large time asymptotic problems because 
the unsteady term vanishes at large time, undesired numerical issues may arise for time-accurate 
computations. One major issue is the error due to stiffness could accumulate after long time integration. 
The boundary conditions must also be treated consistently so that any error or reflection would not produce 
spurious waves that propagate into the computational domain. 

The other issue for most existing finite volume or DG formulations is associated with the Riemann 
solver used to render the flux vector unique at cell interfaces. From the fundamental standpoint, 
formulations for computing flux vectors at the presumed discontinuous cell interface are ad hoc at best. 
Correct (and thus unique) flux vectors at the flow discontinuity can only be approximated accurately for 
any multi-dimensional cell interfaces if and only if the jump solution is known a priori. Any one- 
dimensional based approximation breaks down when the cell interface is not aligned with the discontinuity. 
For this reason, dimensional splitting along with grid adaptation is often needed to capture shocks. Treating 
shock waves at the cell (and flux-integration) interface is apparently departing from fundamental gas 
dynamics. In it, the flux vector is integrated before and after the discontinuity. A control volume contains 
the shock with the boundaries only in the smooth region is used to derive correct jump conditions. While 
such practices are forgotten or neglected by most CFD algorithms, the space-time CESE method is perhaps 
the only method that follows similar concept of flux conservation, i.e., the conservation or control volume 
should be without any discontinuities at the boundaries like most finite-volume family of schemes. 


Following the above-mentioned fundamental concepts, the space-time CESE method is formulated by 


the integral form of the conservation laws with a space-time flux h is defined as: 

h=(F x ,F,F,U) 


( 3 ) 


Based on above-mentioned physical consideration, general conservation laws should be rewritten in the 
following integral form 

j) v h-ds=jSdV (4) 

v 

where the space-time flux vector is integrated over the surface Q of an arbitrary space-time volume V. 
Strictly speaking, the differential form, eq. (1), is only valid for smooth solutions. The surface normal 
vector is defined by s = fl dA , where dA is the area increment on Q, and n is the outward unit normal 


vector. Equation (4) is quite general, and in fact, any conservation law with or without a source vector can 
be cast in this form. Thus, the numerical algorithm devised for eq. (4) can be easily extended to other 
physical problems. Theoretically, the space -time domain can be discretized with time being treated as one 
coordinate. With this unstructured or structured temporal-spatial mesh, a time accurate solution can be 
obtained by integrating eq. (4) using the CESE numerical framework. However, this is not done in most of 
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the CESE implementations due to topological difficulties in generating a four-dimensional mesh for time- 
accurate 3D problem. Alternatively, the time coordinate is organized in a structured, layered fashion with a 
constant increment for either structured or unstructured meshes. 

Given a discretized space-time domain with an unstructured mesh, the dependent variables can be 
assumed to vary according to the following Taylor series expansion truncated to the third order, 

U(x,y,z,t)=U 0 +U t (t-t 0 )+U x (x-x 0 )+U y (y-y 0 )+U z (z-z 0 ) 

+ UJx-x,) 2 /2 + U xy (x-x 0 )(y-y 0 ) + U yy (y-y 0 ) 2 /2 + 

u xz ( x - x o)(z - Zo ) + Uyziy - y 0 )(z ~ z o) + ujz - z 0 ) 2 / 2 + 

ujx - x 0 )(t -t 0 )+ U yt (y - y 0 )(t -t 0 )+ U Zt (z - z 0 )(t - t 0 ) + U„(t -t 0 f / 2 + 

Uxxx (x-x 0 ) 3 /6 + U m (y-y 0 f / 6 + U zzz (z-z 0 f /6 + U m (t-t 0 ) 3 /6 + 
u xxy( x - x o) 2 (y - y 0 ? 1 2 + U^ix - x 0 )(y - y 0 ) 2 / 2 + U xxz (x -x 0 ) 2 (z-z 0 )/ 2+ (5) 

U xzM ~ *o)(* ~z 0 f / 2 + U m (y -y 0 ) 2 (z-z 0 )/2 + U yzz (y - y 0 )(z -z 0 ) 2 / 2 + 
U xxt (x- x o) 2 (. t - t o)^ + U X) , l (x-x 0 )(y-y 0 )(t-t 0 ) + U yyt (y-y 0 ) 2 (t-t 0 )/2 + 
U zzt (z-z 0 ) 2 (t-t 0 )/2 + U yzl (y-y 0 )(z-z 0 )(t-t 0 ) + U xzl (x-x 0 )(z-z 0 )(t-t 0 ) + 

+ U m (x - x 0 )(y - y 0 )(z - z 0 ) + U xtt (x - x 0 )(t - 1 0 ) 2 / 2 + U ytt (y - y 0 )(t -t 0 ) 2 12 + 
U ztt (z-z 0 )(t-t 0 ) 2 / 2 

within each solution element (SE), which roughly coincides with the discretized cell element. The 
coordinate (jc 0 , To, ^o) is a reference point within an SE called solution point in the CESE method. The 
coefficients of the expansion are not equal but analogous to their counterparts in derivatives. These 
coefficients are unknowns to be solved for at each element. Dependent variables vary smoothly within 
each solution element. Jumps are allowed across solution elements. The piecewise smooth assumption 
allows flow discontinuities to be captured while retaining the designed approximations for each element. 
To solve for unknown coefficients, the conservation laws, eq. (4), must be integrated numerically in the 
computational domain. Most high-order unstructured mesh methods carry out numerical integration along 
the element interfaces (which coincides with the cell boundaries) only in the spatial domain. To handle 
jumps across these element interfaces, numerical approximations must be employed to determine unique 
fluxes based on solutions from either side of the interface. In contrast, for the CESE method, integration of 
eq. (4) is not carried out along the boundaries of the solution element, but along a staggered temporal- 
spatial conservation element (CE). The conservation elements are staggered with existing solution 
elements such that all CE interfaces only pass through the smooth region of any solution element. 
Integrating fluxes across a flow discontinuity, similar to what is done to derive Rankine-Hugoniot relations, 
poses no numerical issues at all and no ad hoc Riemann solver is needed. Figure 1 shows the CE 
arrangement around a triangular element. For two-dimensional CESE solvers, one quadrilateral CE is 
constructed for each neighboring element and for 3D tetrahedral elements, the CE for each neighbor is a 
double pyramid. Flux integration using eq. (4) around all surrounding CEs results in surface integrals on 
three different group of faces: i) top faces ( GBAF , GFED , and GDCB ), ii) bottom faces ( GBAF\ GF 
ED , and GDCB), and iii) side faces (all remaining faces). Bottom and side face integrals depend only 
on solutions from previous time level. Top face integration requires dependent variables at the new time 
level. 

If the Taylor series approximations shown in eq. (5) are truncated to first order, the resulting numerical 
method is formally second-order accurate. In the CESE numerical framework, the solution point of each 
SE is located at the geometric center of all surrounding CEs (the union of GBAF , GFED , and GDCB). The 
solution procedure for the second-order accurate schemes involves three steps. 
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5 1 . Integration of eq. (4) can be simplified to algebraic equations containing only unknowns U 0 , at the 
new time level. Due to the choice of solution point, all derivatives ( U x , £7^ and U z ) vanish in these 
equations. Solutions of these algebraic equations can be easily calculated for all elements without 
any matrix inversion. 

52. Coefficients U x , U y and U z for each element are calculated by using generalized finite difference 
formulations for unstructured neighbors. For each neighbor, a set of derivative coefficients is 
calculated by solutions around each edge. 15 A weighted average procedure in favor of smaller 
values (discussed below) is used to determine a set of unique derivative coefficients for each 
element. 

53. Update temporal derivative, U h by using the original governing equations. This completes one time 
iteration. 

Procedure S2 involves spatial derivatives calculations using finite-difference formulations. Unlike 
conventional structured-mesh algorithms, these spatial derivatives are unknowns to be determined at each 
time step in the CESE method. Theoretically, these spatial derivatives can be calculated for each solution 
element by solving flux- conservation equations for each quadrilateral CE. These theoretical base line 
solutions for spatial derivatives form the non-dissipative a-scheme and numerical dissipation can be added 
by modifying them. 13 Several methods can be used to add numerical dissipation including the c-scheme. 13 
Courant number insensitive scheme (CNI), 18 and the edge-based derivative (EBD) method. 15 For highly 
non-uniform or high-aspect ratio meshes, the EBD method has been shown to produce the most robust 
results. 




Figure 1. Conservation elements for a 2D triangular mesh: (a) top view, 
(b) space time CE volume. 


For an unstructured mesh, solutions at three points are needed to determine U x and U y in a two- 
dimensional domain and four points are needed for U x , U y , and U z in a three-dimensional domain. For a 
triangular element with a solution point O and three vertices, V\, V 2 , and U 3; three sets of derivatives based 
on three edges can be computed. As discussed in Ref., 15 these derivatives form an equivalent of the non- 
dissipative scheme. Numerical dissipation can be added by growing the triangles to new vertices , V[ , 
and V[ with a stretching parameter a according to 

v;=o +0^-0) ( 6 ) 

for i = 1, 2, and 3. The parameter, a, has a value greater than one. Larger a implies more numerical 
dissipation. For the edge-based CNI scheme, the o value decreases with the local CFL number. The 
optimum solution accuracy can be obtained with a o value small enough to maintain numerical stability. 
Numerical dissipation control can thus be achieved by decreasing the a value for large cells. For high 
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aspect ratio meshes, the value of a has to be increased significantly to maintain numerical stability. The 
shape of the triangle affects both solution accuracy and numerical stability. For an elongated cell, the edge 
length can vary by orders of magnitude. The shortest side requires a very large a (in the order of 10 3 - 10 5 , 
depending on the aspect ratio) to maintain stability. Stretching by the same ratio for the longest side would 
add excessive amount of numerical dissipation. The combined effect often causes inaccurate viscous 
solution near the wall. Better solution accuracy could be obtained while maintaining numerical stability by 
using a different a for different vertices. Thus, eq. (6) is modified to 

v;~0 + a t {V-0) ( 7 ) 

The stretching parameter, cr i? is determined by 

a,. = max(a min ,a A min /||^-0||+1) (8) 

Additional parameters are introduced to fine tune the control. The global minimum spacing between 
any two vertices is denoted by A min and cr min is a pre-set lower limit of the stretching parameter. The value a 
in eq. (8) is a constant within an element. Its value can be chosen as a global constant or could be a simple 
function of grid sizes (for example, a = CFL cr max , where er max is the maximum stretching desired). For a 
high-aspect ratio mesh, a very large value of a would increase o x significantly only for the short side 
(smaller ||V^ — 0|| values). Numerical experiments indicate that the use of eqs. (7) and (8) allows accurate 

solutions of wall shear stress for a flat-plate boundary layer with an aspect ratio as high as 10 5 . For a 
relatively uniform mesh, o < 2 is sufficient to maintain both accuracy and stability. For a high-aspect ratio 
mesh, on the other hand, a needs to be increased to 10 3 or higher. Nevertheless, such a high value is only 
applied at the vertices with a vey small spacing. 

In general, more numerical dissipation (larger a) is needed in the vicinity of flow discontinuities. A 
shock detection scheme helps locate the flow discontinuities. More numerical dissipation can thus be 
added accordingly. The shock sensor proposed by Ref. 25 is used: 

d p =Vp*(MV) (9) 

where M is the local Mach number and V is the local velocity vector. If S v is greater than a threshold 
value, the element is considered to contain a shock; and proper numerical dissipation can be added 
accordingly. Using non-dimensional values (normalized by freestream values), the value of S v is greater 
than 10 in the vicinity of shocks generally, although this threshold value varies from problem to problem. 
The use of the above shock sensor avoids adding excessive numerical dissipation away from the shock. 
For the high-order CESE method, it can also be used to reduce approximation orders around flow 
discontinuities. 

For triangular and tetrahedral elements, each neighbor provides sufficient information to construct all 
spatial derivatives. Suppose there are n derivatives around an element. The n value varies from 3 for a 
triangle to 12 for a hexahedral element. The weighted averaging procedure by Chang 26 is used. The 
process begin by first calculating 


0,'=V(C) 2 +(C) 2+ (C) 2+£ 

where superscript z indicates the z-th derivative and s = 10" 2 ° is added to avoid zero derivatives. Sort all 6 X 
and define Ci as the sorted array arranged from the largest to the smallest value (i.e., Ci = max (^)). 


n n - 1 T — r n 


By definition, all are positive numbers. The weighting factor for the z-th derivative is then defined as 

W/ = (1+ a/yjCFL Ip, )g, (10) 

where CFL is the local CFL number at each element and a is a global parameter with a range from 0 to 
about 2. Larger a values imply larger numerical dissipation. For acoustic wave simulations, a very small a 
value can be used to increase the phase accuracy. For shock capturing, the a value can be set to 1 or 2. 
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The parameter g t accounts for the geometrical effect of different neighboring cell sizes. The CE volume 
can be assigned to g f such that larger CE carries more weight in the derivative calculations. The unique 
derivatives at each element can be then calculated by 

i = 1 i = 1 

on 

i = 1 i = 1 

i = 1 i = 1 

The above weighted average procedure is used universally for all the computations with or without 
discontinuities. No adjustment or fine-tuning of the weighted averaging is needed to handle bow shocks or 
even flows with very complex shock patterns. This weighted averaging procedure is in a way similar to 
conventional flux limiters used in upwind formulations. Nonetheless, the weighted average procedure in 
the CESE framework does not alter flux vectors directly. The weighted derivative coefficients calculated 
simply modify numerical dissipation near flow discontinuities to maintain numerical stability. 

The second-order CESE methods have been used for steady-state as well as time-accurate computations 
extensively. To further improve numerical efficiency and solution accuracy without significantly increasing 
the total grid counts, the time-accurate local time-stepping and the high-order CESE methods discussed 
below have been investigated in this paper. 


Time-accurate Local Time-stepping CESE Method 

In the initial stage of the CESE method development, an ad hoc local time stepping has been 
implemented for steady state computations. In this attempt, each element is advanced at a different time 
step satisfying the stable CFL constraint (CFL < 1). Flux conservation in time is not preserved. It was 
found that for simple flow geometries such as flat-plate boundary layers, this ad hoc local time-stepping 
procedure does yield rapidly converging and accurate steady-state solutions. However, for flows with 
strong shocks such as blunt body configurations, the calculated shock stand-off distance was wrong 
compared to exact solutions or solutions obtained by constant time step. This is not surprising because 
when flux conservation is not preserved, there is no guarantee for correct steady-state solutions. Chang et 
al., 16 propose a time-accurate local time stepping framework for ID flows by using interface patches. In 
this formulation, flux conservation in both space and time is ensured by book keeping solutions at multiple 
time steps. It has been shown that solution accuracy is preserved even with a large step size jump across 
the cell interface. In this investigation, we follow similar procedure to implement time-accurate local time- 
stepping (TALTS) for multi-dimensional flows. 

The fundamental concept of TALTS is based on space-time flux conservation. Conservation in the 
spatial domain alone does not guarantee flux conservation in time, a key ingredient for temporal accuracy 
with disparities in time steps. The CESE method is constructed based on flux conservation in the temporal- 
spatial domain; and, thus is fully equipped to implement TALTS. The following algorithms are taken to 
implement TALTS in the Navier-Stokes code: 

LI. Sorting time steps: local time steps are first calculated based on a prescribed CFL number 
(normally between 0.5 and 1). Assume the smallest time step found is At, an array of time steps, 
f(k) = 2 k At is constructed. The local time step for each element is adjusted to the nearest /(£) value. 
After the adjustment, time steps vary as a power of 2 throughout the entire domain. For a relatively 
uniform mesh, the maximum k value needed is no more than 3 or 4. For viscous flow 
computations, the maximum k value could be as large as 10 or more. A maximum k value of 10 
amounts to three orders of magnitude variation in time steps (or element sizes). RANS calculations 
could require a maximum k of 15 or more. 
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L2. Determining solution levels: base on the time steps of neighbors, each element must determine how 
many solution levels must be stored for flux integration. Normally 2 or 4 levels of solutions need to 
be stored. If the cell size jumps abruptly, more levels are required. For a smooth mesh, no more 
than two levels are needed. 

L3. Integrating flux with patches: for cells with smaller time-step neighbors, flux integration on the 
side faces must be done with patches. Solutions stored with multiple time levels are used to 
construct flux vectors on the divided side-face patches. Summing flux vectors for all the patches 
completes the side-face integration to ensure flux conservation in time and space. 

L4. Timing iterations: a clock to monitor physical time is added for time marching. An element will 
march to the next time level only when the local cell time after the iteration matches the current 
physical time. Caution must also be taken to avoid advancing time before the neighbors catch up. 

L5. Communicating multiple solution levels: for parallel computations, the intermediate solutions must 
be communicated to neighboring domains for accurate flux integrations. 

The baseline algorithm in the ez4d code has been modified slightly for viscous flow computations. 15 The 
original strong flux conservation formulation in conjunction with a modified edge-base derivatives 
approach is used to effectively compute viscous flow with highly stretched unstructured meshes inside the 
boundary layer. Solution accuracy has been carefully validated for time-accurate inviscid and viscous 
flows. The ez4d code is developed in the C++ language entirely by following modern software engineering 
procedures. Object-oriented and generic programming using templates have been employed extensively in 
the software framework for ez4d. Sound software architecture ensures easy maintenance and agile future 
development. Both message passing interface (MPI) and multi-thread capabilities are integrated in the 
software framework to take advantage of multi-core, distributed memory cluster environments. Currently, 
the ez4d code can solve Euler or laminar Navier-Stokes equations for triangular, quadrilateral, tetrahedral, 
or hexahedral unstructured meshes. 

The 4th-Order Accurate CESE method for 3D Navier-Stokes Equations 

The implementation of the 4th-order CESE method for 3D Navier-Stokes equations follows the 
theoretical framework introduced by Chang for ID scalar equations 18 . Traditionally, high-order methods 
are implemented by increasing the order of the truncated Taylor series representation of a solution element. 
These high order polynomials are used for the integration of the governing equations. While the treatment 
of integrations varies from finite volume, spectral element, spectral differences to DG, the unknown high- 
order coefficients are obtained by solving a matrix equation for all flow dependent variables. Increasing the 
order of the polynomials is equivalent to sub-dividing an element into smaller sub-elements. The 
additional polynomial coefficients, or quadrature points, decrease the size of the domain of dependence. 
As a result, the allowable time step decreases significantly as the order of polynomials increases. For 
example, the maximum CFL number is only around 0.1 for a typical 4th-order DG method. 19 The reduced 
time step size partially offsets the increasing accuracy of the solutions with fewer elements. To avoid such 
a stringent time step requirement, an implicit integration technique can be used. However, such implicit 
method is difficult to implement for unstructured meshes and only works for steady-state problems. The 
4th-order accurate CESE method introduced by Chang, on the other hand, preserves the maximum CFL 
number around 1 for explicit formulations. The formulation of this 4th-order scheme for Navier-Stokes 
equations will be discussed below with emphasizes on multi-dimensional and general flux formulations. 

The fundamental building block for preserving the CFL limit in the high-order CESE method is the 
introduction of additional equations for higher derivatives. For a 4th-order scheme, equations for the 
second derivatives of the dependent variables can be obtained by differentiating the governing equations 
twice. Six second-derivative equations for Navier-Stokes equations can be derived as follows. 
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For simple conservation laws such as the linear wave equation, the second-derivative equations take an 
identical form as the original governing equation. To solve these second-derivative equations, additional 
space-time flux vectors are introduced: 

h = (F F F U ) 

,l xx Jtxx ’ i xxy ’ ± xxz ’ xx / 


h =(F F F U ) 

,l xy j ixy’ x xyy’ x xyz’ xy > 

h =(F F F U ) 

,l yy xyy ’ 1 yyy’ 1 yyz’ ^ yy > 

h =(F F F U ) 

h =(F F F U ) 

h =(F F F U ) 

n zz ^ xzz ’ JZZ’ 1 zzz’ U zz' 


(13) 


For compressible Navier-Stokes equations, a total of six conservation laws listed below need to be 
integrated at each time level in addition to eq. (4). 

<£ h rr • ds =0 

J a(V) xx 


(£ h xv -ds= 0 

Ja(V) xy 

(£ h -ds= 0 

Ja(v) yy 

(£ h Y7 -ds= 0 

J a(V) xz 

(£ h V7 -ds= 0 

Ja(V) yz 

<£ h„-ds= 0 

J a(V) zz 


(14) 


The Mathematica® 27 commercial software is used in this paper to derive all high- derivative governing 
equations and to generate corresponding C++ code to be included in the ez4d code. For 3D flows, six 
double derivatives, U xx , U xy . U m U yz , U xz , and U zz , need to be solved in addition to U. The solution 
procedure for the 4th-order scheme can be summarized as follows: 

HI. Solving second derivatives: using eq.(14) and the procedure SI used for the second-order method, 
all six second derivatives are calculated for each element. 
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H2. Computing third spatial derivatives: using the computed second derivatives and the procedure S2 
for the second-order method, third spatial derivatives are computed for each element. 

H3. Integrating high-order flux equations: using the computed second and third derivatives, solution U 0 
can be computed by integrating eq.(4) by using eq.(5). Note here that with the selection of the 
solution point being the geometric center of the combined CEs, the first derivatives still vanish in 
the top faces flux integrals. Second and third derivatives on the top faces are already determined in 
HI and H2. This results in a simple and easy way to solve algebraic equations for C/ 0 . 

H4. Calculating first spatial derivatives: first derivatives are calculated using a procedure similar to S2 
but with high-order finite difference formulations. The same weighted averaging as in the second- 
order scheme is applied (eq. (11)). 

H5. Updating temporal derivatives: using all spatial derivative coefficients computed, all derivatives 
with temporal differentials are calculated using the governing equations. This completes the 
iteration. 


The high-order scheme requires some additional equations and geometric information to be computed. 
Temporal partial derivatives (calculated in H5) are obtained by differentiating the governing equations, for 
example, 


u tt 


dx dy dz 


(15) 


All corresponding derivatives on the right-hand side can be obtained by analytical expressions by using 
the derivative coefficients defined in eq. (5). Again, all these temporal partial derivatives are derived by 
using the Mathematica® software to avoid errors. In procedures H3 and H4 above, the integration of the 
top and bottom faces requires high-order moments associated with integrations of the second and third 
order terms in eq. (5). These high-order moments are pre-computed and stored for each CE with respect to 
the CE center before the iteration starts. There are total of seven moments required for 2D calculations and 
16 for 3D geometries. 

In the procedure H4, eq. (4) contains integrations of flux vectors over the side faces. These flux vectors 
are functions of the high-order polynomials defined in eq. (5). To maintain high-order accuracy, Jacobians 
of these flux vectors with respect to U, and all derivatives must be derived and used for flux integration 
over space and time if the flux vectors are evaluated as a function of polynomials directly. Deriving these 
Jacobian terms for 3D Navier-Stokes equations with variable viscosity involves tedious book keeping 
processes and is thus error-prone. In this paper, an alternative procedure is used. All high-order flux 
integrations on the side faces are carried out using quadrature integrals. 28 Flux vectors are calculated at 
quadrature points within the side face by using the high-order solution element already computed. The flux 
integral is then computed using high-order quadrature formulae. No Jacobian terms are needed in this 
approach. The nine-point Gaussian quadrature is used for 2D elements and a combined triangular 4-point 
quadrature in space and 3 -point in time is used for 3D elements. These quadrature points are designed to 
carry out integrals over the quadrilateral and triangular prism side faces, respectively. 


III. Results and Discussion 

The time-accurate local time-stepping and the 4th order CESE methods discussed above are aimed at 
providing higher-fidelity numerical solutions for unsteady flow computations. Several test cases are 
investigated in this paper to assess the effectiveness of both approaches. In the next two subsections, these 
test cases are described in detail. 

Time- Accurate Local Times-stepping for Unsteady Flows 

The TALTS method is designed to handle large grid size disparities in a discretized domain. One 
frequently encountered problem is to simulate acoustic waves propagating through a viscous mesh. 
Computing acoustic waves using such a non-uniform mesh is quite challenging. Figure 2 shows a 
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rectangular domain with a clustered mesh near the bottom boundary. A constant time step computation 
results in significant phase error due to the large CFL number disparities. The wave amplitude is also 
incorrectly damped. In contrast, the TALTS computation maintains a rather uniform CFL number (around 
0.8) and thus results in much improved accuracy in both amplitude and phase throughout the domain. 

Flux conservation plays a crucial role in shock capturing schemes. It is well known that failing to 
conserve the flux in the numerical computations could result in a wrong solution. An ad hoc local time- 
stepping method uses different time steps at different cells and flux conservation in time is violated. 
Conventional schemes use local time stepping for steady-state flow computations in which the residual 
norm is driven to a small value at convergence. For the CESE method, such an approach may work 
differently for steady state computations. The violation of flux conservation in time persists for long time 
integration. As a result, the converged solutions only guarantee a balance between temporal and spatial 
fluxes and may not preserve the desired flux conservation at large time. One example investigated here is a 
Mach 6 flow over a cylinder computed with a typical wall-clustered viscous mesh. The mesh has a 
maximum aspect ratio of about 10 3 and tenth-power local time steps (k = 10) were used. Figure 3 shows 
the computed Mach number contours compared against inviscid shock-fitted solutions 29 using two different 
methods of time marching. The ad hoc local time-stepping results clearly over-predict the shock stand-off 
distance. The TALTS results on the other hand correctly compute the shock locations. The Mach number 
contours in the inviscid region also agree well with shock-fitting Euler solutions. The ad hoc local time- 
stepping method, despite its drawbacks in shock capturing, handles smooth steady-state solutions with good 
accuracy. It also has the characteristics of quick convergence to large-time asymptotic solutions. The best 
strategy for steady-state computations is to first run ad hoc local time-stepping to quickly settle with a 
nearly correct solution. The TALTS method is then used to converge the solution to the correct large-time 
asymptotes. 


Near- Wall Stretched Mesh 




0.99992 0.99995 0.99998 1.00001 1.00004 


Time-accurate local time-stepping 



x 


Figure 2. Acoustic waves propagating through a domain with clustered mesh, showing density 
contours at a time instant. 
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Figure 3. Mach 6 viscous flow over a cylinder with a wall-clustered viscous mesh, Mach number 
contours compared inviscid solutions obtained by shock-fitting method: 29 (a) ad hoc local time- 
stepping (b) TALTS. 

In more practical applications, a large mesh-size disparity is unavoidable. Flow physics dictates the 
required mesh size. The ideal mesh should have a grid distribution that is commensurate with the physical 
length scale. The development of the TALTS scheme is aimed at coping with a computational domain with 
various physical length scales. In the next test case, a cone with a counterflowing jet at the leading edge 
submerged in a Mach 1.6 free stream is computed. The shear layer and jet instability waves alter the 
leading edge shock. To accurately compute the flowfield for sonic boom applications, both the near- field 
instability waves and far- field shocks need to be captured with proper resolutions. To reduce the total grid 
counts, the mesh has to transition from the near field fine grid to the farfield coarse grid like that shown in 
Fig. 4. In the constant time step approach, the finest mesh near the jet exit determines the time step size 
used in the simulations by using the constraint of CFL < 1 . In the farfield, the same time step implies a 
very small local CFL number. Numerical dissipation increases significantly when the CFL number 
decreases for the CESE method. The computed normalized pressure disturbance contours (based on the 
freestream pressure value) shown in Fig. 5(a) exhibit a twist near the grid transition region. Further 
downstream, the shock appears to be diffused due to excessive numerical dissipation. On the contrary, the 
computed solutions using TALTS displayed in Fig. 5(b) show consistent shock resolutions across both 
regions. In this particular example, the grid size varies several orders of magnitude. The TALTS 
simulations use £=11, i.e., the ratio between the largest and the smallest time step is 2,048. 

Another application of the TALTS method is to compute steady-state supersonic flow over a large blunt 
body as shown in Fig. 6. To resolve the viscous boundary layer and the strong bow-shock simultaneously, 
the mesh size varies by several orders of magnitude. The TALTS scheme can significantly improve the 
accuracy both in shock-capturing and resolving boundary- layer profiles. In this particular example, ninth- 
power time levels ( k = 9) for the 32 million tetrahedral elements were used. Parallel efficiency for the 
TALTS schemes needs some careful consideration. The mesh in this example is decomposed into 120 
blocks running in parallel using 120 cores on NASA Langley Research Center’s K-cluster with Sandy 
Bridge Intel CPUs. The average physical time required is about 1.38 seconds per time iteration. Similar 
computations using constant time step took about 4.2 seconds physical time per time iteration. This 
amounts to about a factor of 3 saving in computational time by using the TALTS method. The speed-up in 
computational efficiency is reasonable. However, load balancing is suboptimal because many blocks only 
contain relatively large cells. These blocks would spend significant amount of time waiting for other 
blocks with small time steps to catch up. To further improve parallel efficiency, the domain decomposition 
must be done with due consideration of equal grid size distribution. A “chunk” style grid with each chuck 
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including roughly equal amount of fine and coarse elements can significantly improve overall throughput 
for large-scale computations. 



(a) (b) 


Figure 4. Mesh transition from fine near-field to coarse farfield used for sonic boom applications: (a) 
mesh interface between two zones (pressure contours levels the same as Fig. 5(a) below (b) enlarged 
view of the triangular mesh. 



Figure 5. Computed pressure contours for a Mach 1.6 flow over a cone with a counterflowing jet 
located at x = 1: (a) constant time step (b) TALTS. 
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Figure 6. Computed Mach number contours at two planes for a Mach 6 flow over a blunt body with 
a 28° angle of attack, computed using TALTS. 


The 4th-order CESE Method for Scalar Equations 

To verify the order of accuracy for the 4th-order CESE numerical framework, the following scalar wave 
equation is used for extensive tests. 


du ,du du du x 

— + a ( — + — + — ) = 0 (16) 

dt dx dy dz 


The corresponding flux vector is h = ( au , au , au , u). Without loss of generality, the parameter value for 
a is set to be 1. Due to the simplicity of this equation form, the high-order equations (eq. (13)) for the 
above scalar equation take the same form as the original equation. Equation (16) allows a variety of exact 
solutions to be constructed for testing purposes. A polynomial solution 


u=m, e=t~ (i7) 

is used for validation. For simplicity, the function is assumed to satisfy f{6) = 6 n , where n is a positive 
integer. A 20x20^20 cubic domain is discretized using four different meshes: 10 3 , 20 3 , 40 3 , and 80 3 . The 
resulting structured meshes are further sliced into tetrahedrons. The above exact solutions are imposed as 
initial and boundary conditions. Numerical dissipation parameters are fixed at a = 0 and a = 1.5. 
Theoretically, when the polynomial order satisfies n < 4, the 4th-order scheme should give exact solutions. 
And indeed, the computed solutions retained exact solutions with an error norm hovering around 10" n -10" 15 
for a large number of iterations. When ^ = 4, the error norm versus number of grid counts in each direction 
is shown in Fig. 7. As expected, the error norms converge at 4th-, 3rd-, and 2nd-order rate for U 0 , U X9 U xx , 
and Uxxx, respectively. The error norms for y - derivatives converge at an identical rate. 
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Figure 7. L2 error norms of U 0 and 
derivatives versus 1/N (N is the number of 

points in each direction) for / ( 6 ) = 6 4 (Lines 
correspond to theoretical 4th-, 3rd-, and 
2nd-order convergence). 


Figure 8. L2 error norms of U 0 and 
derivatives versus 1/N (N is the number of 

points in each direction) for f(6) = Asin(0) 
(Lines correspond to theoretical 4th-, 3rd-, 
and 2nd-order convergence). 


The second test is done for an exact wave solution with the following form: 

f(6)=A w sin(/9), d=cot-k x x-k y y-k z z 

The relation between the frequency and wave numbers is (O — k x — k y — k z =0 To fully test the multi- 

dimensional formulations, the wave angle is 45° with respect to the x-axis. The computed L2 error norms 
for dependent variables and their derivatives are shown in Fig. 8. Similar 4th-order convergence is verified 
for the wave solutions. 

Acoustic wave propagation through a square 2D domain is used to validate the 4th-order accurate 
algorithms for 2D Euler equations. Isotropic triangular meshes as shown in Fig. 9(a) are used for accuracy 
studies. Four different mesh levels with 462, 1,880, 7,378, and 29,276 elements are used for computations. 
These meshes correspond to an average cell size decreases of factors 2, 4, and 8, respectively. Exact 
solutions of linear acoustic waves are imposed both as initial and boundary conditions. The wave 
propagates along the diagonal direction as shown by the contour plots in Fig. 9(b). In each case, the L2 
error norm is computed after 20 periods of waves have passed through the domain. A CFL number of 0.8 
is used for all four levels of meshes. Numerical dissipation parameters are set to a = 0 and <7 = 2. Lower 
numerical dissipation by decreasing a can increase the accuracy and lower the error norm. However, for 
order of accuracy studies, a constant value of 2 is used. 

Exact solutions of acoustic waves are obtained by assuming linear perturbations. Error norms for all 
four meshes are first calculated using linear Euler equations with a density acoustic wave amplitude of 10." 4 
Figure 10(a) shows the expected slightly better than the 4th-, 3rd-, and 2nd-order convergence for 
dependent variables and their derivatives, respectively. The same set of computations is repeated using the 
full Euler equations, and the results are shown in Fig. 10(b). Surprisingly, except for the third derivatives, 
the error norms of all other variables and derivatives level off around 10." 8 A literature search confirmed 
that most, if not all, acoustic wave validation cases for high-order methods had been performed for the 
linearized Euler equations because the exact acoustic wave solutions are only valid for linear equations and 
solutions calculated by full Euler equations have nonlinear effects. However, a simple estimate can be 
made. For a wave amplitude of 10," 4 quadratic nonlinear terms modify the linear solutions at a value 
around 10" 8 . Thus, the error norm can never be smaller than 10" 8 . To verify this, the acoustic wave 
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amplitude is changed to 10" 6 and the resulting norm is shown in Fig. 11(a). The L2 error norms now 
converge at the correct rate for dependent variables and all derivatives. Nonlinear terms only affect the 10" 
12 level, lower than the error of the finest mesh error norm. If the amplitude is further dropped to 1 O' 8 , the 
nonlinear terms affect only at the machine accuracy level at 10" 16 . Consequently, the error norms converge 
at the expected rate without any problems as shown in Fig. 1 1(b). 



(a) (b) 

Figure 9. Acoustic waves propagating diagonally through a 20 x 20 square domain using isotropic 
triangular meshes (a) the coarsest mesh (b) velocity contours indicating wave direction. 




(a) (b) 

Figure 10. L2 error norm versus number of elements, N, for acoustic waves propagating through a 
square domain with a density amplitude of 10" 4 : (a) linear Euler solutions (b) full Euler solutions. 
Lines indicate theoretical 4th-order convergence. 
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Figure 11. L2 error norm versus number of elements N for acoustic waves propagating through a 
square domain with different density amplitude: (a) full Euler solutions with amplitude = 10" 6 (b) full 
Euler solutions with amplitude = 10." 8 Lines indicate theoretical 4th-, 3rd-, and 2nd-order 
convergence. 

The test cases established by Shu and Osher 22 for ID acoustic and shock wave interaction has been used 
extensively for testing high-order algorithms. The shock impinges on an acoustic wave and generates 
short-scale oscillations along with shocklets upstream. Numerical computations must capture the shock 
and preserve the penetrated acoustic wave without excessive damping. The shock detection method is used 
to switch to second-order schemes in the vicinity of shocks (the dissipation parameter a is increased to 5 at 
the shock). The computational results using 1,600 points by Shu and Osher has been treated as the grid- 
converged solution in many code validation cases. The 4th-order CESE method is also tested with the 
same number of points. Figure 12(a) compares the computed density distribution at one time instant. The 
good agreement is evident. Both the shorter scale acoustic waves near the shock and N-shaped shocklets 
further upstream are properly resolved by the CESE method. The agreement in phase (locations of peaks 
and valleys) is also quite remarkable. However, the CESE results predict slightly larger (about 5-10%) near 
shock acoustic waves amplitudes. Figure 12(b) shows a similar calculation but with a lower weighted- 
average dissipation parameter, a = 0.5. The computed penetrated acoustic wave amplitudes increase even 
further. Additional computations using 3,200 points indicate similar phenomena. Results calculated by the 
4th-order CESE method do not appear to converge at the consensus grid-converged solution for this 
benchmark problem. The consensus grid-converged solution represents only one discretized solution using 
a relatively larger amount of numerical dissipation for the CESE method. Reducing this numerical 
dissipation results in an increase of the penetrated acoustic wave amplitudes. Due to the lack of 
experimental data, the true converged solution is not yet determined. Further computations are needed to 
shed more light on this benchmark problem. 
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Figure 12. Comparison of computed density distribution using 4th-order CESE with Shu and 
Osher 22 by using a o value of 1.5 and different numerical dissipation: (a) a = 1 (b) a = 0.5. 

As a more difficult test case involving complex geometry and shock patterns, the Mach 1 .4 flow over a 
multi-gutter wall, previously investigated experimentally by Suzuki et al. 24 and computationally by Wang 
et al., 25 is tested with the high-order CESE method. A total of about 125,000 triangular elements were used 
in the computations. Figure 13 compares results obtained by using both second- and 4th-order CESE 
methods using the same mesh with the experiments. The shock sensor defined in eq. (9) is used to detect 
shocks. For both computations, slightly more numerical dissipation (<r = 2) is used in the vicinity of shocks 
and a uniform dissipation parameter (cr = 1.1) is used elsewhere. In the 4th-order computations, the 
discretization is also switched to second-order in the vicinity of shocks. Keeping high-order derivatives 
across the shock is numerically unnecessary. Switching to second-order help resolve the shocks better. As 
can be seen from Fig. 13, both computations give quite satisfactory shock resolutions. The 4th-order 
solutions provide more details around the gutter wall and near the complex reflected shocklets, additional 
short-scale waves reflecting from the wall is also visible. This example demonstrates the use of high-order 
method in the region where small-scale physics are present. Numerical computations are also underway to 
detect the minimum grid requirement for the high-order CESE method. 
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(C) 


Figure 13. Normal density gradient contours for a Mach 1.41 flow over a multi-gutter wall: (a) 
results with second-order CESE (b) results with 4th-order CESE (c) experimental Schlieren 
photography by Suzuki et al. 24 at a slightly different time. 


As a preliminary step for validating the developed Navier-Stokes solvers, a Mach 3 flow over a flat 
plate is computed by using a triangular mesh sliced from a structured 41x21 quadrilateral grid. The inflow 
and top boundaries are fixed at a Mach number of 3, the bottom wall boundary is set to be at adiabatic 
conditions, and the downstream boundary uses non-reflecting conditions. The computed velocity contours 
overlay by the mesh used is shown in Fig. 14(a). The velocity profiles at a Reynolds number of 10 5 agree 
well with the compressible boundary solutions, as shown in Fig. 14(b). With only about 10 mesh points 
within the boundary layer, the agreement is satisfactory for a 4th-order Navier-Stokes solver. More 
detailed validations including order of accuracy are under investigation and will be the subject of a separate 
paper. 


0.05 | 



(a) 



(b) 


Figure 14. Mach 3 viscous flow over an adiabatic flat-plate boundary layer computed by high-oder 
CESE: (a) velocity contours and grid distribution (b) velocity profiles at Re = 10 5 compared with 
compressible boundary-layers solutions. 
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IV. Concluding Remarks 


As an research effort toward developing next-generation, high-fidelity time-accurate flow solvers for 
complex geometries, the space-time CESE numerical framework has been extended for two new numerical 
algorithms. The TALTS method allows different time steps to be employed at different cells while 
maintaining flux conservation over the entire space-time domain. Traditional explicit time-accurate solvers 
suffer from stringent time step restrictions caused by the finest mesh near either the viscous wall or minute 
geometrical details. A great deal of computational time is wasted at large cells during the time 
advancement. For schemes like CESE that are sensitive to the local CFL number, the solution accuracy 
also degrades at large cells due to extremely small time steps. The TALTS scheme eliminates these 
numerical issues for problems with a large grid size disparity. Numerical results indicate that in addition to 
better computational time compared to the conventional constant time-step simulations, the solution 
accuracy is greatly enhanced for the CESE framework. This new capability has the potential to help future 
computations for problems with multiple scales and physics such as transitional or turbulent flow 
simulations. Load balancing could be a potential issue for the TALTS schemes if the decomposed domains 
do not contain evenly distributed gird size distributions. However, more careful domain decompositions 
could alleviate this issue. 

The CESE flow numerical framework has also been extended to allow 4th-order accurate Navier-Stokes 
computations. Solutions obtained for all 2D test cases show excellent accuracy with designed order of 
convergence. The stable CFL number limit remains around one for the 4th-order accurate CESE method. 
This feature eliminates the need to decrease the time step size substantially as the order of accuracy 
increases. Scalar wave, Euler, and Navier-Stokes equations for scalar waves, acoustic waves, and 
boundary-layer flows are validated for simple test cases. The flow variables, first-, send-, and third- 
derivatives have L2 error norms that converge at 4th-, 3rd-, and 2nd-order, respectively. The high-order 
framework can thus be used in the regions where derivatives are more important. The test case for a 
supersonic flow over a multi-gutter wall shows improved resolution away from the shock by using the 
high-order CESE method. The combination of the TALTS schemes and the high-order CESE framework 
could potentially provide a new computational paradigm for future large-scale time-accurate simulations. 
More validations for 3D cases and the combination of TALTS with high-order methods are underway, and 
results will be reported in a future paper. 
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